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Tlie Fermi-Pasta-Ulam (FPU) model, which was proposed 50 years ago to examine thermalization 
in non-metallic solids and develop "experimental" techniques for studying nonlinear problems, con- 
tinues to yield a wealth of results in the theory and applications of nonlinear Hamiltonian systems 
with many degrees of freedom. Inspired by the studies of this seminal model, solitary-wave dynamics 
in lattice dynamical systems have proven vitally important in a diverse range of physical problems — 
including energy relaxation in solids, denaturation of the DNA double strand, self-trapping of light 
in arrays of optical waveguides, and Bose-Einstein condensates (BECs) in optical lattices. BECS, 
in particular, due to their widely ranging and easily manipulated dynamical apparatuses — with one 
to three spatial dimensions, positive-to-negative tuning of the nonlinearity, one to multiple com- 
ponents, and numerous experimentally accessible external trapping potentials — provide one of the 
most fertile grounds for the analysis of solitary waves and their interactions. In this paper, we 
review recent research on BECs in the presence of deep periodic potentials, which can be reduced 
to nonlinear chains in appropriate circumstances. These reductions, in turn, exhibit many of the 
remarkable nonlinear structures (including solitons, intrinsic localized modes, and vortices) that lie 
at the heart of the nonlinear science research seeded by the FPU paradigm. 

PACS numbers; 05.45.Yv, 03.75.Lm, 03.75.Nt, 05.30.Jp 



The Fermi-Pasta-Ulam (FPU) model was for- 
mulated in 1954 in an attempt to explain heat 
conduction in non-metallic lattices and develop 
"experimental" (computational) methods for re- 
search on nonlinear dynamical systems fl*]. Fur- 
ther studies of this problem ten years later led 
to the first analytical description of solitons (us- 
ing the Korteweg-de Vries equation, which is a 
continuum approximation of the discrete FPU 
system), which have since become one of the 
most fundamental paradigms of nonlinear sci- 
ence. These nonlinear waves occur ubiquitously 
in rather diverse physical situations ranging from 
water waves to plasmas, optical fibers, super- 
conductors (long Josephson junctions), quantum 
field theories, and more. Over the past several 
years, the study of solitons and coherent struc- 
tures in Bose-Einstein condensates (BECs) has 
come to the forefront of experimental and theo- 
retical effbrts in soft condensed matter physics, 
drawing the attention of atomic and nonlinear 
physicists alike. Observed experimentally for the 
first time in 1995 in vapors of sodium and ru- 
bidium [2, 3], a BEC — a macroscopic cloud of co- 
herent quantum matter — is attained when (10^ to 



10^) atoms, confined in magnetic traps, are opti- 
cally and evaporatively cooled to a fraction of a 
microkelvin. The macroscopic behavior of BECs 
near zero temperature is modeled very well by 
the Gross-Pitaevskii equation (a time-dependent 
nonlinear Schrodinger equation with an external 
potential), which admits a wide range of coher- 
ent structure solutions. Especially attractive is 
that experimentalists can now engineer a wide 
variety of external trapping potentials (of either 
magnetic or optical origin) confining the conden- 
sate. As a key example, we focus on BECs loaded 
into deep, spatially periodic optical potentials, ef- 
fectively splitting the condensate into a chain of 
linearly-interacting, nonlinear droplets, the dy- 
namics of which is accurately characterized by 
nonlinear lattice models. This paper highlights 
some of the quasi-discrete nonlinear dynamical 
structures in BECs reminiscent of the discover- 
ies that originated from the FPU model. 
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I. INTRODUCTION 

One of the most important nonlinear problems, whose 
origin dates back to the early 20th century, concerns the 
conduction of heat in dielectric crystals. As early as 1914, 
Peter Debye suggested that the finite thermal conductiv- 
ity of such lattices is due to the nonlinear interactions 
among lattice vibrations (i.e., phonon-phonon scattering) 
. To understand the process of thermalization — which 
refers to how and to what extent energy is transported 
from coherent modes and macroscopic scales to internal, 
microscopic ones — and to develop computational tech- 
niques for studying nonlinear dynamical systems, Fermi, 
Pasta, and Ulam (FPU) posed the following question in 
1954: How long does it take for long-wavelength oscil- 
lations to transfer their energy into an equilibrium dis- 
tribution in a one-dimensional string of nonlincarly in- 
teracting particles? This question has since spawned a 
diverse array of activities attempting to answer it and 
fomented a strong impetus to research in topics such as 
soliton theory, discrete lattice dynamics, and KAM the- 
ory. Furthermore, these fronts remain active research 

topics HIE Hi, 13. 

Before the FPU work, it was commonly assumed that 
high-dimensional Hamiltonian systems behave ergodi- 
cally in the sense that a smooth initial energy distribu- 
tion should quickly relax until it is ultimately distributed 
evenly among all of the system's modes (that is, thermal- 
ization should occur). To explicitly verify this fundamen- 
tal hypothesis of statistical physics, FPU constructed a 
one-dimensional dynamical lattice, with N identical par- 
ticles which interact according to an anharmonic repul- 
sive force. 

Running numerical simulations on the computers avail- 
able in the early fifties, FPU observed that the lattice did 
not relax to thermal equilibrium, contrary to everybody's 
expectations P, H El m m. An especially striking 
observation was a beating effect, in the form of a near- 
recurrence of the initial long-scale configuration, which 
reappeared after a large number of oscillations involving 
short-scale modes. In this manner, more than 97% of 
the energy returned to the initial mode. Moreover, this 
finding was robust with respect to variations in the total 
number of particles and particular choice of the (power- 
law) anharmonicity in the interaction between them. 

Motivated by this study, Zabusky and Kruskal con- 
sidered a continuum version of the model, showing that 
the dynamics of small-amplitude, long-wavelength per- 
turbations obeys [on the timescale ^ (wavelength)"^] the 
Korteweg-de Vries (KdV) equation. They subsequently 
introduced the concept of solitons (solitary waves) in 
terms of the KdV equation 0, IT^ ll^ . The explana- 
tion for the lack of thermalization is that the energy 
gets concentrated in robust coherent structures (the soli- 
tons), which interact elastically and thus do not trans- 
fer their energy into linear lattice modes (phonon waves, 
also referred to as "radiation") Q. The KdV equation 
thereby became the first example in the celebrated class 



of nonlinear partial differential equations (PDEs) that 
are integrable by means of the inverse scattering trans- 
form. It was later followed by numerous other impor- 
tant PDEs, such as the nonlinear Schrodinger equation, 
the sine-Gordon equation, higher dimensional examples 
(including the Kadomtsev-Petviashvili and the Davey- 
Stcwartson equations) and multi-component examples 
(including the Manakov equation) J^j . 

Since then, the study of solitons, and more general co- 
herent structures, has become one of the paradigms of 
nonlinear science |12L ITH . Il7| | . Such dynamical behavior 
occurs in a wide variety of physical systems — including 
(to name just a few examples) nonlinear optics, fluid me- 
chanics, plasma physics, and quantum field theory. Over 
the past several years, the impact of solitary-wave dy- 
namics has been especially significant in the study of 
Bose-Einstein condensates (BECs) In this short 

review, we focus on this application and, in particular, 
its description in physically appropriate cases in terms of 
dynamical lattices. 

The rest of this paper is organized as follows. In Sec- 
tionm we define the FPU problem and briefly survey its 
mathematical properties. In Section IIIII we provide an 
introduction to Bose-Einstein condensates and their soli- 
tary wave solutions. We consider, in particular, BECs 
loaded into optical lattices (OLs), and use a Wannier- 
function expansion to derive a dynamical lattice model 
describing this system. In appropriate limits, this leads 
to a (generally, multiple-component) discrete nonlinear 
Schrodinger (DNLS) equation 20). In Section lTVI we ex- 
amine the regime in which a deep, spatially periodic OL 
potential effectively fragments the BEC into a chain of 
weakly interacting droplets. The resulting model, which 
consists of a Toda lattice with on-site potentials, pro- 
duces self-localized modes that may be construed as soli- 
tons of the underlying BEC. Section W\ concludes the pa- 
per, discussing a number of future directions. 



II. THE FPU PROBLEM 

The FPU model consists of a chain of particles con- 
nected by nonlinear springs. The constitutive law of the 
model, i.e., the relation between the interaction force F 
and the distance y between adjacent atoms in the chain 
(with only nearest-neighbor interactions postulated), was 
taken to be F{y) = —[y + G{y)]. FPU considered three 
different functional forms of G{y): quadratic, cubic, and 
piecewise linear (lOi] . 

In the case of a cubic force law ^iv) = ^{v + 
3/3y'^), where f3 is an effective anharmonicity coefficient, 
one may write 



(1) 



vj = + L(%+i-yj) 

+ ivj - yj-i f + ivj+i - yj){yj - v^-i)] | > 

which is supplemented with fixed boundary conditions, 
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yo = Vn = 0. FPU 0, 01 used the initial condition 
yj{0) = sinO'Tr/TV) and yj{0) = 0. This form CJ of 
the FPU chain is obtained by discretizing the continu- 
ous nonlinear string, 



d^y 



1 + 3/3 



dy_ 

dx 



(2) 



with the following approximations to the continuum 
derivatives: 



dx 
dx^ 



vj+i - yj yj - yj-i 
h ' h 



yj+i - "^yj + yj~i 



(3) 



Here, yj is the displacement of the j-th particle from its 
equilibrium position, h = L/N is the normalized spac- 
ing (the distance between the particles), L is the string's 
length, and N is the number of particles, which FPU took 
to be 16, 32, or 64 in their calculations. (Using equation 
the onset of resonance overlaps was studied in the 
FPU problem in the first ever application of Chirikov's 
overlap criterion 

With an appropriate scaling {t ht, y ^ hy^Ji/ (3), 
the FPU chain can also be written 



% = - {yo ^vj-i) 

+ \ [ivj+i - yjf - ivj - yj-if] 

Using the continuum field variable u |22j. 

1 rv^ 



yt_ 

'2h 



(4) 



(5) 



Eq. |(5J) yields, to lowest order in h, the modified 
Korteweg-de Vries (mKdV) equation (with r = h^t/24:, 
^ = X — ht), 



12v?u^ + u^^^ — , 



(6) 



which is further reduced to the KdV equation proper via 
the Miura transformation 0, 0| . 

One can also derive the KdV equation directly from an 
FPU chain with a quadratic anharmonicity in the inter- 
particle interaction [l3,|ll|ll, F{y) = -{y + ay^). In 
the latter case, the discrete FPU model takes the form 

^, = ^^±l^^|±^(l + c.^^±i^). (7) 

To study the near-recurrence phenomenon in Eq. |(7J, 
Zabusky and Kruskal H El derived its continuum limit 
{h — >0,Nh — > 1), 

ytt = yxx + eyxyxx + -^h^yxxxx + 0{eh'^,h^) , (8) 
where e = 2ah. 



Unidirectional asymptotic solutions to Eq. (jSJ are con- 
structed with (l^i 

2/^0(e,T), ^ = x-t, T = ^et, (9) 
where the function (j) obeys the equation 

+ kki + + O {h\h^e-^) = , (10) 



for small 5^ = h'^e ^/12. Finally, with u 
obtains the KdV equation from Eq. ifTUI) 23]: 

Ur + UU(^ + S'^U^^^ — . 



one 



(11) 



Eq. Hll() has solitary-wave solutions of the form 

u{x, t) = 2K2sech^ [k{x - Ak^I - xq)] , (12) 

with constants k and xq. [Note that a more rigorous 
derivation of the KdV equation from the FPU chain (as 
a fixed point of a renormalization process) has recently 
been developed Q.] 

Although the solitary-pulse solution H12|l has been well- 
known since the original paper by Korteweg and de Vries, 
it was the paper by Zabusky and Kruskal '13] that re- 
vealed the particle-like behavior of the pulses in numer- 
ical simulations. (The term "soliton" was coined in that 
paper to describe them.) Since then, solitons have be- 
come ubiquitous, as their study has yielded vital insights 
into numerous physical problems 0, 0, 0, 0, 0, |3 ■ 
In the next section, we discuss their importance to Bose- 
Einstein condensation 0, 0| . 

It is remarkable that even today, 50 years from the orig- 
inal derivation jlQi] , FPU chains are themselves still stud- 
ied as a means of understanding a variety of nonlinear 
phenomena. Recent studies focus not only on the model's 
solitary- wave solutions and their stability ^^^Sn •25i | , 
but also on its thermodynamic properties and connec- 
tions with the Fourier law of heat conductivity [31 . its 
dynamical systems/invariant manifold aspects l2 6|. l27l . 
and its connections with weak turbulence theory |28l | . 

To conclude this section, we remark that the most nat- 
ural discrete model that has been derived from the NLS 
(or GP) equation for a soliton train is the Toda lattice 
|22(| (see also the details discussed below). However, the 
leading-order nonlinear truncation of the latter lattice 
equation once again yields the FPU model. Conversely, 
one can approximate the FPU chain by the NLS equation 
in the high-frequency limit 0, [s^] ■ The validity of this 
approximation varies with time due to energy exchange 
between modes. 



III. BOSE-EINSTEIN CONDENSATION 

At low temperatures, bosonic particles in a dilute gas 
can reside in the same quantum (ground) state, forming 
a Bose-Einstein condensate (BEC) El El IM • Sev- 
enty years after they were first predicted theoretically. 
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BECs were finally observed ejcperimentally in 1995 in va- 
pors of rubidium and sodium [^|^ . In these experiments, 
atoms were loaded into magnetic traps and evaporatively 
cooled to temperatures on the order of a fraction of a mi- 
crokelvin. To record the properties of the BEC, the con- 
fining trap was then switched off, and the expanding gas 
was optically imaged |l9j |. A sharp peak in the velocity 
distribution was observed below a critical temperature, 
indicating that condensation had occurred. 

Under experimental conditions, BECs are inhomoge- 
neous, so condensation can be observed in both momen- 
tum and coordinate space. The number of condensed 
atoms Af ranges from several thousand (or less) to several 
million (or more). The magnetic traps confining BECs 
are usually well-approximated by harmonic potentials. 
There are two characteristic length scales. One is the 
harmonic oscillator length, a^o — ■^/^/(rnwho) (which is, 
typically, on the order of a few microns), where m is 
the atomic mass and Who — {'-^x^y^zY^^ is the geomet- 
ric mean of the trapping frequencies. The s econd scale 
is the mean healing length, x = l/y^STrjaln, where n is 
the mean density of the atoms, and a, the (two-body) 
s-wave scatteringlength, is determined by collisions be- 
tween atoms [i^ [l^ 1^ ,3^ . Interactions between atoms 
are repulsive when a > and attractive when a < 0. The 
length scales in BECs should be compared to those in 
condensed media like superfluid helium, in which the ef- 
fects of inhomogeneity occur on a microscopic scale fixed 
by the interatomic distance 19]. 

With two-body collisions described in the mean-field 
approximation, a dilute Bose-Einstein gas is very accu- 
rately modeled by the cubic nonlinear Schrodinger equa- 
tion (NLS) with an external potential [i.e., by the so- 
called Gross-Pitaevskii (GP) equation]. An important 
case is that of cigar-shaped BECs, which are tightly con- 
fined in two transverse directions (with the radius on the 
order of the healing length) and quasi-free in the longitu- 
dinal dimension 3.5, 36. 37, 38, 39]. In this regime, one 
employs the ID limit of the 3D mean-field theory, gen- 
erated by averaging in the transverse plane (rather than 
the ID mean-field theory per se, which would be appro- 
priate were the transverse dimension on the order of the 
atomic size [H HI IH H ) . 

The original GP equation, describing the BEC near 
zero temperature, is Isl Isl lioHilf 

= + + ^(^^)) ^^^^ 

where ^I' = ^'(r, t) is the condensate wave function nor- 
malized to the number of atoms A/", V{r) is the external 
potential, and the effective interaction constant is 
50 = [iTih^a/m][l + 0{C)l where C = \/\WW is the 
dilute-gas parameter. The resulting normalized form of 
the ID equation is pi IMIstLIsI jsol] 

i^t = -\-^xx + g\ip?i' + y{x)i', (14) 

where ■0 and V are, respectively, the rescaled ID wave 



function (a result of averaging in the transverse di- 
rections) and external potential. The rescaled self- 
interaction parameter g is tunable (even its sign), because 
the scattering length a can be adjusted using magnetic 
fields in the vicinity of a Feshbach resonance (42, ||4J] • 

A. BECs in Optical Lattices and Superlattices 

BECs can be loaded into optical lattices (or superlat- 
tices, which are small-scale lattices subjected to a long- 
scale periodic modulation), which are created experimen- 
tally as interference patterns of counter-propagating laser 
beams HHHIillillilliSllSIIll. Over the past sev- 
eral years, a vast research literature has developed con- 
cerningBECs in suchpotentials [H HI H El El El 
l58l Isg . l60l l6lL |62, |63, IM, |65j , as they are of considerable 
interest both experimentally and theoretically. Among 
other phenomena, they have been used to study Joseph- 
son effects 52], squeezed state s ISq . Landau-Zener tun- 
neling and Bloch oscillations '53], controllable conden- 
sate splitting [sH . and the transition between superflu- 
idity and Mott insulation at both the classical [b^, IH] 
and quantum E^ levels. Moreover, with each lattice site 
occupied by one alkali atom in its ground state, BECs in 
optical lattices show promise as a register in a quantum 
computer |6^l70l]. 

With the periodic potential V{x) — V{x + L), one may 
examine stationary solutions to (|14() in the form 

il;{x,t) ^ R{x)eyi^{i[e{x) ~ ^Jit]) , (15) 

where /i, the BECs chemical potential, is determined 
by the number of atoms in the BEC; it is positive for 
repulsive BECs and can assume either sign for attrac- 
tive BECs. Using the relation dd/dx — c/B? ("angular 
momentum" conservation), one derives a parametrically 
forced Duffing-oscillator equation for the amplitude func- 
tion HiilHllzilzlIzI, 

R" + f^R - gR^ ~ v{x)R = , (16) 

where R" = d^R/dx^. 

Equation H16|l admits both localized and spatially ex- 
tended solutions. Supplemented with appropriate bound- 
ary conditions, it yields both bright and dark solitons, 
which correspond, respectively, to localized humps on the 
zero background, and localized dips in a finite-density 
background. These states are similar to the bright and 
dark solitons in nonlinear optics; they are stable, respec- 
tively, in attractive and repulsive ID BECs |74Ll75l ]. 

When V{x) is spatially periodic, the bright solitons 
resemble gap solitons, which are supported by Bragg 
gratings in nonlinear optical systems. In BECs, they 
have been observed in two situations: (1) the small- 
amplitude limit, with the value of /x close to forbid- 
den zones ("gaps") of the underlying linear Schrodinger 
equation with a periodic potential (7^ Ft^ : and (2) in 
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the tight-binding approximation (discussed below), for 
which the continuous NLS equation can be replaced by 
its discrete counterpart — the so-called discrete nonlinear 
Schrodinger (DNLS) equation "ts]. In the latter context, 
the strongly localized solutions arc known as intrinsic 
localized modes (ILMs) or discrete breathers. Spatially 
extended wave functions with periodic or quasi-periodic 
R{x), which may be either resonant or non-resonant with 
respect to the periodic potential V{x), are known as mod- 
ulated amplitude waves and have been shown to be sta- 
ble (against arbitrary small perturbations) in some cases 

ill mill in III- 



B. Lattice Dynamics 

In the presence of a strong optical lattice, the GP 
equation H14|l can be reduced to the DNLS equation 
[23,13123 • To justify this approximation, the wave func- 
tion is expanded in terms of a set of Wannier functions 
localized near the minima of the potential wells. 

The eigenvalue problem associated with the linear part 
of Eq. ^ is 



(17) 



where ipk,a can be expressed in terms of Floquet-Bloch 
functions, (pk,a — e''^^Ufc,Q,(x), with Uk,a{x) — Uk,a{x+L) 
and a indexing the energy bands, so that Ea (k) = Ea (fc+ 
[27r/L]) [s^. The energy is represented using Fourier 
series, 



E^{k) 



iknL 



where the asterisk denotes complex conjugation and 



(18) 



2tt 



r/L 



Ea,{k)e 



-iknL 



dk. 



(19) 



-■k/L 



Although the Floquet-Bloch functions provide a com- 
plete orthonormal basis, it is more convenient to utilize 
Wannier functions (also indexed by a), 

Wo.{x - nL) - I ' (^fe,„(x)e^™'=^dfc , (20) 

V 27r 

which are centered about x — nL (n G Z). The Wannier 
functions constitute a complete orthonormal basis with 
respect to both n and a. One can also guarantee that 
the Wannier functions are real by conveniently choosing 
phases of the Floquet-Bloch functions. 

Given the orthonormality of the Wannier function ba- 
sis, any solution of H14|l can be expanded in the form 



ip{x,t) 



(21) 



with coefficients satisfying a DNLS equation with long- 
range interactions, 

-dCn.a 



dt 



rii 



^ y / , / , ^ni,Qi"-n2,«2'^n3:a3 "^'0,01,02,03 

01,02,03 ni,n2,Ti3 

(22) 

as was illustrated in Q^. In (j^ . 

/C30 
U'n,o'U^rn^ajlUn2,02^n3,03<^2; . (23) 
-oo 

Because the Wannier functions are real, the integral in 
(|23|1 is symmetric with respect to all permutations of both 
(a, ai, Q!2, as) and (n, ni, n2, ns). 

Although Eqs. (|22|l are intractable as written, sev- 
eral important special cases can be studied ^20]. The 
nearest-neighbor coupling approximation is valid when 
\^\.a\ ^ |w„^Q.| for n > 1. More generally, one can assume 
that the Fourier coefficients in H18|) decay rapidly beyond 
a finite number of harmonics. This simplifies the lin- 
ear term in (|22|l . Additionally, because Wn,a is localized 
about its center at a: = nL, it is sometimes reasonable to 
assume that the coefficients satisfying n = ni = n2 = n^ 
dominate W^'^^llallal^ so that the others may be ne- 
glected. In the nearest-neighbor regime, this implies that 



dr 



dt 



9 ^ ^ W^O,Oi,02,03C„ Q,jC„^Q.2C„.ct3 
01,02,03 



(24) 



where Wa,aua2,a3 = W^Q,o;"o"o3 independent of n 
ic tight-binding model, 

'^0,oCn,o + '<^l,o(Cn-l,Q + Cn+l.o) 



This leads to the tight-binding model, 
.dc 



dt 



9W0. 



(25) 



in the single-band approximation. Typically, this approx- 
imation is valid if the height of the barrier between po- 
tential wells is large and if the wells are well-separated. 
While this intuition may be generally true, the Wannier 
function reduction provides a systematic tool that can 
establish the validity of the approximation on a case by 
case basis (by determining the overlap coefficients); see, 
e.g., |23| for specific examples. Including next-nearest- 
neighbor coupling in this regime allows one to study in- 
teractions between intrasite and intersite nonlinearities. 

In more general situations in which single-band 
descriptions are inadequate because of resonant, 
nonlinearity-induced interband interactions, one can sim- 
plify by assuming that Coi^a = wi,Q,(e) ~ 0{e) as e ^ 
0. Applying the phase shift Cn,a{t) — exp[ia)o,oi]cno(i), 
so that the nonlinear terms consist of oscillatory expo- 
nents, and supposing that c„^q,(0) is sufficiently small 
and may be ignored, we conclude that, on the time 
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scale ^ 1/e, the exponents oscillate rapidly except when 
a = a2, ot\ — a^i or a — a-i^ a\ = a-i- Through time- 
averaging, one obtains 



dt 



'^l,a{Cn-l,a + Cn+l,a) 
"1 



(26) 



where Wa^ai = Wa.ai,a,ai dcscribcs the interband in- 
teractions. Equation H26|l is a vector DNLS with cross- 
phase-modulation nonlinear coupling. An example of the 
implementation of this method is illustrated in Fig. ^ 
More generally, the advantage of the approach of Ref. 
[20| is that, given the explicit form of the potential, the 
relevant coefficients can be computed and the appropri- 
ate reduced (single-band or multiple-band) model can be 
derived to the desired level of approximation. 







FIG. 1: Comparison of the lattice reconstructed solution in 
the tight-binding (dashed line) and the 3-band (dash-dotted 
line) approximation with the exact solution (solid line). The 
comparison is performed for V{x) = — 5cos(2a;) and different 
chemical potentials: fi = 1.5 (left panels) and /i — —1.5 (right 
panels). The bottom (semi-log) panels show the result of 
dynamical time-evolution of the tight-binding approximation 
with the dotted line. As time evolves (we depict snapshots 
at f « 50), the latter can be seen to approach the shape 
of the exact solution (in the left panel, these plots cannot 
be distinguished) and to match its asymptotic form, possibly 
shedding small wakes of low-amplitude wave radiation in the 
process (see, for example, the bottom right panel). 

At this stage, one can study the ILMs of Eqs. H25|l or 
(|26|l and use the Wannier function expansion to recon- 
struct the solution of the original GP equation |23| . This 
approach has been successfully used in a variety of ap- 
plications and for different localized states present in the 
lattice — including bright, dark, and discrete-gap solitons, 
as well as breathers |63, U3 uR IHI ■ Other phenomena, 
such as discrete modulational instabilities, have also been 
studied [t^ • More specifically, one of the most successful 
implementations of discrete NLS equations (and variants 



thereof) in this context included the quantitative predic- 
tion that its modulational stability analysis determined 
the threshold of a dynamical instability of the condensate 
(the so-called classical superfiuid-insulator transition of 
|67|). These predictions were subsequently verified quan- 
titatively by the experimental measurements of Ref. |68| . 



IV. SOLITON-SOLITON TAIL-MEDIATED 
INTERACTIONS AND THE TODA LATTICE 

Recent advances in trapping techniques allow the gen- 
eration of bright solitons and chains of bright solitons in 
effectively ID attractive condensates 82, 83, 84, 85] . In 
this section, we consider the collective motion of a chain 
of bright solitons. We focus, in particular, on the dynam- 
ics of attractive BECs trapped in a deep optical lattice 
(OL) that renders a ID attractive condensate into a chain 
of interacting solitons (see Fig. [2). 

Consider a BEC loaded into an OL potential pro- 
duced by the interferencepattern of multiple counter- 
propagating laser beams dlt AS, 49j 50, ^ij. 
In principle, it is possible to design various optical trap 
profiles (essentially at will) by appropriately superimpos- 
ing interference patterns. For the purposes of this expo- 
sition, we adopt an OL profile, with a tunable inter-well 
separation, given by the Jacobi elliptic-sine function. 



V{x) = Vbsn2(x;fc), 



(27) 



where Vq is the strength of the OL. The elliptic modulus 
k allows one to tune the separation between consecutive 
wells, r = ^oy+i - ^0;j = '^K{k), where ^o;i = 2j K{k) 
is the position of the j-th well and K{k) is the complete 
elliptic integral of the first kind. 
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FIG. 2: A quasi-lD condensate (solid line) in a deep peri- 
odic optical- lattice potential (dashed line). The condensate 
is effectively described as a chain of coupled solitons whose 
positions follow a Toda lattice with on-site potentials [Eq. 
I|31|l ]. Using the oscillating ansatz l|82|l . where the j-th soli- 
ton is forced to oscillate with amplitude Aj, one can further 
reduce the dynamics to a second-order recurrence relationship 
between neighboring amplitudes [Eq. Ij^H^ ]. 



The stability properties of BECs in the o ptic al lattice 
potential f?7| have been recently studied |62l F^ . W^ . 
Here, we are interested in the case of large separation 
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between the wells {k ~ 1), when the BEC is effectively 
fragmented into a chain of nearly identical solitons with 
tail-mediated interactions, subject to the action of an 
effective on-site potential due to the OL. For a large set 
of parameter values, the system can be reduced to a Toda 
lattice as was illustrated in Refs. [H The 

reduction involves two steps. 

First, the interaction between the j-th soliton [at posi- 
tion {t)] and the j-th well is approximated, using vari- 
ational techniques js^ or methods based on conserved 
quantities js^, by 

= -Kff(^.-eo;,), (28) 

which describes a particle in the effective potential Vcff 
felt by the soliton. For well-separated troughs (fc — 1), 
the effective potential may be approximated by |88j| 

/ 224 \ 

K'ff(0 « ^Vo (^aiC - j^^Vo + a,ej , (29) 

where v is the average amplitude (height) of the soliton, 
ai = 8/15 and 03 — —16/63. 

Second, one treats the interaction of consecutive soli- 
tons in the absence of the OL. This interaction is well- 
studied in the context of optical solitons |2^ |8^ [o^, loil 
For identical, well-separated solitons (with 
the phase difference tt between adjacent solitons), it is 
approximated by the Toda-lattice equation for the soli- 
ton positions, 

= TL{Cj-i,£,j,£,j+i) 

(30) 

Finally, after combining (|30|) with the on-site potential 
dynamics (|28|l , we reduce the dynamics of a weakly cou- 
pled BEC in the deep OL to the Toda lattice with on-site 
potentials, 

= Tl{^j-1,Cj , 6-+l) - Kff (6- - Co;,) . (31) 

The Toda lattice H3()ll admits exact traveling-soliton 
solutions However, the effective potential in ll^ 

breaks the lattice's translational invariance and accom- 
modates the existence of ILMs (breathers). To describe 
such localized oscillations, we consider small vibrations 
about the equilibrium state {^j = Co;j)- 

=Cou cos(wt), (32) 

where uj is the common oscillation frequency for all soli- 
tons, and the j-th soliton vibrates with an amplitude 
Aj about its equilibrium position ^q.j (see Fig.|2l. Sub- 
stituting the ansatz H32|l into Eq. (|31|l and discarding 
higher-order modes yields a recurrence relationship be- 
tween consecutive amplitudes, 

An+l ^ {a + b Al) An ~ Ar^^l , (33) 



where a = 2 - cu"^ + iaid, b = Sa^d, d = e'^7(16z^'*), r 
is the separation between adjacent troughs, and ai and 
03 are the coefficients from the expansion of the effective 
potential H29|l . Note that the method just described is 
applicable to any OL profile that reduces the dynamics 
of a single soliton to that of a particle inside an effective 
potential given by a cubic polynomial in ^. 

By defining consecutive oscillation amplitudes as j/„ = 
An and x„ = An-i, the recurrence relationship (|33|l can 
be cast as a 2D map, 

f Xn+l = Un 

(34) 

I Vn+i = {a + byi)yn- Xn 

It is crucial to note here that the iteration index n in 
(|34|l corresponds to the spatial lattice-site index for the 
soliton chain. Thus, forward (backward) iteration of the 
2D map H34|) amounts to a right (left) shift of the spatial 
position in the soliton chain. As our goal is to find spa- 
tially localized states of the BEC chain (i.e., ILMs), we 
are interested in finding solutions of Eq. (|33|l for which 
A„ ^ as n ^ ±00 and > (see Fig. That 
is, {xn,yn) (0,0) as 71 — + ±00. These solutions corre- 
spond to homoclinic connections of the origin. 

A typical intersection for the stable (solid curve) and 
unstable (dashed curve) manifolds of the 2D map H34I) is 
depicted in the top- left panel of Fig. |21 The intersection 
points Pi between the stable and unstable manifolds then 
generate a localized configuration for the recurrence re- 
lationship H33() , as depicted in the top-right panel of Fig. 
O When this localized configuration is inserted into the 
original CP equation H14I) . one obtains a spatially local- 
ized, multi-soliton state (depicted in the bottom panel 
of Fig. O by shifting each soliton from its equilibrium 
position by the prescribed amount f87,88l. 

It is important to note that generic localized ini- 
tial configurations do not give rise to long-lived, self- 
sustained ILMs. Nonetheless, the construction described 
above is quite efficient in producing approximate initial 
configurations that generate clean and robust localized 
states, such as the one displayed in Fig. |21 The struc- 
tural and dynamical stability of these localized states is 
quite interesting. For example, the structural stability 
of the homoclinic tangle of the 2D map guarantees the 
existence of the ILM solution in the original model H14|l . 
despite the employment of various approximations in the 
former's derivation. On the other hand, the dynamical 
stability of ILMs permits their observation even in the 
presence of strong perturbations (such as noise). Indeed, 
numerical experiments show that ILMs prevail even in 
the presence of large perturbations to the initial config- 
uration or strong numerical noise, as discussed in Refs. 

mil. 

Finally, it is worth mentioning that the dynamical 
reduction to the 2D map H34() can also be used to 
generate — in addition to the localized states discussed 
above — a wide range of spatiotemporal structures by fol- 
lowing the map's fixed points, periodic orbits, quasi- 
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FIG. 3: Homoclinic connection of the origin (top-left panel) 
[from Ref. |8^1. giving rise to a spatially localized profile (top- 
right panel) . Bottom: the localized state in the original BEG 
model 1141 generated by the prescribed amplitude configura- 
tion. The shaded base, depicting d\u\/dt, highlights the areas 
in which the atomic density \u(x, t)\ varies the most. Observe 
that the solution decays as one moves away from the center 
(x = 0). (As the solution is symmetric with respect to a:: = 0, 
only its right-hand side is shown.) 



periodic orbits, and even chaotic orbits. Further, the 
techniques described in this section can also be applied 
to chains of bright solitons in which the deep OL is re- 
placed by an array of focused laser beams or impurities 
that tends to pin the solitons and serve as a local effective 
attractive potential [see Eq. if^ ] . 



ence. Specifically, using Wannier-function expansions, we 
reduced the dynamics of the Gross-Pitaevskii (GP) equa- 
tion to a discrete nonlinear Schrodinger (DNLS) equa- 
tion, which was subsequently used to identify localized 
solutions (as well as a number of other interesting dy- 
namical instabilities, such as those discussed in Ref. 79[) 
and reconstruct the nonlinear waves of the original GP 
equation. Treating a BEG trapped in a strong OL as a 
collection of interacting solitons, we then used perturba- 
tive and/or variational techniques to reduce the original 
GP equation to a Toda-lattice equation with an effective 
on-site potential. This nonlinear lattice supports robust, 
exponentially localized oscillations in BEG soliton chains. 

The crucial ingredient, amenable to the dynamical sys- 
tems perspective described in this manuscript, is the GP 
equation's nonlinear self-interaction, which is induced as 
a mean-field representation of two-particle interactions 
and is responsible for a plethora of dynamical behaviors. 
The most remarkable of these, solitary waves (solitons), 
lie at the heart of the paradigm that originated from the 
seminal work of Fermi, Pasta, and Ulam. 

There remain, moreover, a multitude of exciting di- 
rections for future research. For example, the study of 
nonlinear lattice models in higher dimensions has yielded 
a number of interesting generalizations and more^exotic 
solitary waves such as multi-dimensional solitons _95] and 
discrete vortices [o^ . A systematic analysis of their 
existence properties, dynamical characteristics (mobility 
and structural stability), thermodynamic properties, and 
physical relevance involves a number of interesting and 
subtle questions that may keep nonlinear scientists like us 
busy until we're ready to celebrate the 100th anniversary 
of the publication of the FPU problem. 



V. CONCLUSIONS 

In this work, we have surveyed recent research on lat- 
tice dynamics of Bose-Einstein condensates (BECs) in 
optical lattice (OL) potentials. We discussed, in partic- 
ular, the dynamics of BEGs in a deep OL that naturally 
discretizes the original spatially continuous system into 
a collection of discrete interacting nodes (wells). The en- 
suing nonlinear lattice features the trademark localized 
solutions — solitons and intrinsic localized modes — that 
have rendered the Fermi-Pasta-Ulam problem one of last 
century's most important contributions to nonlinear sci- 
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